/**
   Application.cxx
   \file Definition of class Application
 
   \author Thomas Burton
   \date 12/6/11
   \copyright 2011 BNL
*/

#include "tppmc/Application.h"

#include <cmath>

#include <TROOT.h>
#include <TCollection.h>
#include <TNamed.h>

#include <TClonesArray.h>
#include <TDirectory.h>
#include <TMath.h>
#include <TPythia6.h>
#include <TTree.h>

#include "tppmc/DistributionQaFigures.h"
#include "tppmc/Distributions.h"
#include "tppmc/Enumerations.h"
#include "tppmc/Event.h"
#include "tppmc/FragmentationFunction.h"
#include "tppmc/PartonDistribution.h"

namespace tppmc {

  Application::Application()
    : mTree(NULL)
    , mEvent(NULL) {
    Splash(std::cout);
  }

  /** \todo add try...catch in case of bad configuration */
  Application::Application(Configuration& config)
    : mConfig(config)
    , mTree(NULL)
    , mEvent(NULL) {
    Splash(std::cout);
    Initialise(config);
  }

  Application::~Application() {
    TIter next(gROOT->GetListOfFunctions());
    // Most pointers are held by smart pointers.
    // TTree is owned by the output file and is deleted when
    // the file closes, so we don't delete it.
    if(mEvent) {
      delete mEvent;
      mEvent = NULL;
    } // if
  }

  void Application::Splash(std::ostream& os) const {
    os << std::endl <<
      "                 TPPMC" << std::endl <<
      "Transverse Polarized Proton Monte Carlo" << std::endl <<
      "              Version " << smVersionMajor << "." << smVersionMinor <<
      std::endl << std::endl <<
      "Monte Carlo Generator for asymmetry studies\n" <<
      "   in polarized proton-proton collisions\n\n" <<
      "      http://code.google.com/p/tppmc/" <<
      std::endl;
  }

  bool Application::Initialise(Configuration& config) {
    if(not config.Verify()) {
      return false;
    } // if
      // Suppress printing of messages about pages added to postscript,
      // new canvases etc
    gErrorIgnoreLevel = 5000;
      
    mConfig = config;
    std::cout << "\nHere is your run configuration:\n" << std::endl;
    config.Speak(std::cout);
      
    mRandom.reset(new TRandom3(0));
      
    mPdfs.reset(new Distributions(config));
    std::cout << "The distribution functions to be used are:" << std::endl;
    mPdfs->Speak();
      
    mTreeFile.reset(new TFile(config.outfile.c_str(), "RECREATE"));
    if(mTree) {
      delete mTree;
    } // if
    mTree = new TTree("events", "TPPMC event tree");
    if(mEvent) {
      delete mEvent;
    } // if
    mEvent = new Event(0., config.EtaMin(), config.EtaMax(),
		       config.MinHadronPt(),config.tparticleout);
    mTree->Branch("event", &mEvent, 32000, 99);
      
    if(config.GenerateQaFigures()) {
      std::cout << "Generating qa figures to " << config.QaFiguresFileName() << std::endl;
      TDirectory* dir = gDirectory; // Remember initial directory
         
      mQaFile.reset(new TFile(config.QaFiguresFileName().c_str(),
			      "RECREATE"));
         
      tppmc::DistributionQaFigures figures(NULL);
      figures.Fill(*mPdfs);
      figures.Write();
         
      // Restore the initial directory as the current directory
      gDirectory = dir;
      if(dir) dir->cd();
      std::cout << "Finished QA figures" << std::endl;
    } // if
      // Creating a TStopwatch also starts it, so we do this last -
      // we immediately go into the event loop after Initialise() finishes.
    mTimer.reset(new TStopwatch);
    return true;
  }

  bool Application::ProcessEvent() {
    /**
       \todo Protect status message from printing multiple times for a given event
       count if there are cuts on events (eta, pt).
    */
    //timer and counter will appear further down instead
    /*      if(mTree->GetEntries() % 100 == 0 and mTree->GetEntries() not_eq 0) {
    // Print the time since the timer began
    std::cout << "Generated " << mTree->GetEntries() << " events in "
    << mTimer->RealTime() // Stops the timer
    << " seconds" << std::endl;
    // Restart the timer
    mTimer->Start(true);
    }*/ // if
    // Generate a vertex for this event 
    if((mConfig.machine_type not_eq FNAL_E704_PROTON) and
       (mConfig.machine_type not_eq FNAL_E704_PBAR)) {
      double zvtx(0.);
      do {
	zvtx = mRandom->Gaus(0.0, 23.);
      } while(std::fabs(zvtx) > 30.);
    } // if
      // Generate a pythia event and get the particles
    TPythia6::Instance()->GenerateEvent();
    // Rotate the event by a phi angle
    // (pythia generates scattered partons at fixed phi)
    mEvent->SetRotation(mRandom->Uniform(-TMath::Pi(), TMath::Pi()));
    int hadronsize =    mEvent->Update(*mPdfs);
    // Fill the tree if we have hadrons.
    // If we are not fragmentating (function = 0) always fill, because
    // we never have hadrons but we want to save the parton info.
    //    if(mEvent->Hadrons().size() > 0 or mConfig.FragmentationFunction() == 0 or mEvent->Particles().size()>0 ) {
    if(hadronsize > 0 or mConfig.FragmentationFunction() == 0 ) {
      mTree->Fill();
    } // if
    return true;
  }

  bool Application::Close() {
    TPythia6 *pythia = TPythia6::Instance();
    pythia->Pystat(1);
    Double_t xsection = pythia->GetPARI(1);
    UInt_t ntrials = pythia->GetMSTI(5);
    TTree *tree = new TTree("Xsection","Pythia cross section");
    tree->Branch("xsection", &xsection, "X/D");
    tree->Branch("ntrials" , &ntrials , "X/i");
    tree->Fill();
    mPdfs->Unpolarised()->Write("unpol");
    mPdfs->Transversity()->Write("transversity");
    mPdfs->Sivers()->Write("sivers");
    mPdfs->Fragmentation()->Write("frag");
    return mTreeFile->Write() > 0;
  }

  bool Application::EventLoop(int n) {
    if(n < 1) {
      n = mConfig.NEvents();
    } // if
    int ntried=0;
    while(mTree->GetEntries() < n) {
      ProcessEvent();
      if(ntried % 100 == 0 and mTree->GetEntries() not_eq 0) {
	// Print the time since the timer began
	std::cout << "Generated " << ntried<< " for " << mTree->GetEntries() << " selected events in "
		  << mTimer->RealTime() // Stops the timer
		  << " seconds" << std::endl;
	// Restart the timer
	mTimer->Start(true);
      } // if


      ntried++;

    } // while
    std::cout << "tried " << ntried << " events for " << n << " selected events" << std::endl;
    return mTree->GetEntries() == n;
  }
} // namespace tppmc
